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CENTRAL  DIFFERENCE  TIME  INTEGRATOR 


Abstract 


This  paper  is  a  user's  guide  for  the 
stand-alone  explicit  direct  time  integration 
package  STINT/CD  for  structural  dynamics  an¬ 
alysis.  STINT/CD  uses  an  automatic  variable 
time  increment  central  difference  method. 
The  purpose,  function,  limitations,  and 
usage  of  the  package  are  described.  A  FOR¬ 
TRAN  listing  of  STINT/CD  is  given  along  with 
a  sample  problem  which  illustrates  its  usage 
and  performance. 


1.  Introduction 

Direct  time  integration  of  the  discrete  equations  of  motion 
governing  linear  and  nonlinear  structural  dynamics  is  a  fre¬ 
quently  used  solution  method  for  transient  response  analysis. 
It  would  appear  that  the  central  difference  integrator  is  the 
most  commonly  used  explicit  method;  for  example,  it  is  used  in 
the  HONDO  (1]  and  STAGS  [21  computer  codes.  A  difficulty  with 
the  central  difference  integrator  has  been  in  the  selection  of 
the  time  increment  to  maintain  stability  and  desired  accuracy 
over  a  range  of  structural  behavior  and  loading  conditions.  A 
central  difference  integration  method  [3,4]  has  recently  been 
developed  to  overcome  this  difficulty.  This  central  difference 
method  (3,4]  automatically  selects  the  time  increment  to  main¬ 
tain  stability  and  to  achieve  user  requested  accuracy.  The  time 
increment  is  selected  so  that  user  specified  sampling  rates  with 
respect  to  the  dominant  and  maximum  "apparent”  response  frequen- 
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cies  are  maintained.  In  addition  the  central  difference  method 
[3,4]  is  implemented  as  a  stand-alone  package  (STINT/CD)  so  that 
it  is  easily  interfaced  with  existing  structural  analyzers  (fin¬ 
ite  -element  and  -difference  computer  codes) . 

This  paper  is  a  user's  guide  for  STINT/CD  (Stand-alone  Time 
INTegrator/Central  Difference)  and  includes  a  FORTRAN  listing 
(Appendix  A).  The  user's  guide  presents:  1)  a  description  of 
this  automatic  variable  time  increment  central  difference  meth¬ 
od,  2)  a  description  of  the  user  written  subroutines  (examples 
are  listed  in  Appendix  B)  required  to  interface  with  a  host 
structural  analyzer,  3)  a  sample  problem  (embedded  in  the  exam¬ 
ple  user  written  subroutines)  which  illustrates  many  of  the  fea- 

t 

tures  of  this  method,  and  4)  some  suggestions  for  usage  of 
STINT/CD.  The  subroutines  are  written  in  FORTRAN  IV,  utilize 
in-corj  storage  after  an  initial  data  transfer  from  mass  sto¬ 
rage,  and  are  operational  with  minor  changes  on  the  DEC 
VAX-11/780,  UNIVAC  1100,  and  CDC  6000/7000  computer  operating 
systems. 


2.  Method  Description 

The  purpose  and  function  of  this  central  difference  method 
are  described  briefly;  to  the  extent  that  a  user  can  apply  it 
to  a  specific  problem.  The  theory  and  implementation  underlying 
method  are  contained  in  references  (3,4);  the  user  is  en- 


the 
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couraged  to  read  these  references  before  using  the  method.  In 
the  function  description  the  error  measures  are  reviewed  because 
there  are  additions  since  publication  of  references  (3,411.  This 

**  L 

£  section  concludes  with  a  description  of  the  limitations  of  the 
method. 

Purpose 

The  time  integrator  package  STINT/CD,  which  is  comprised  of 
the  subroutines  STINTC,  CENDIF,  ACLTN,  ERRORE,  STEPSZ,  ROTATE, 
and  VI PDA  given  in  Appendix  A,  solves  the  discrete  equation  of 
motion 

M  iin  +  D  +  f  (u1 H)  =  f(tn)  (1) 

“  ***  “■  “S 

governing  structural  dynamics,  where  u  is  the  computational  dis¬ 
placement  vector,  the  superscript  n  indicates  the  vector  value 
at  the  discrete  time  tn,  the  superscript  dot  •  denotes  temporal 
differentiation,  M  is  a  diagonal  mass  matrix,  D  is  a  damping  ma¬ 
trix,  fg(un)  is  the  set  of  internal  forces  due  to  the  stiffness 
that  oppose  the  structural  deformation  (stiffness  forces) ,  and 
f(tn)  is  the  applied  load.  For  a  linear  problem  the  stiffness 
force  fg(u)  becomes  K  u,  where  K  is  the  linear  stiffness  matrix. 
The  damping  term  is  treated  as  _f^«  D  u  ,  or  D  itself  may  be  used 
if  it  is  diagonal.  Therefore  all  quantities  in  eqn.  (1)  are 
computational  vectors;  i.e.  one  dimensional  arrays  of  length 
equal  to  the  number  of  degrees  of  freedom. 
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Function 

The  subroutine  STINTC  provides  an  internal  control  inter¬ 
face  between  the  host  structural  analyzer  and  the  subroutine 
CENDIF  which  is  the  "main"  subroutine  for  the  variable  time  in¬ 
crement  central  difference  integrator.  The  input  to  STINTC  is 
supplied  by  the  user  through  the  user  written  subroutine  DRIVER 
which  is  described  in  section  3.  In  addition  the  diagonal  M  ma¬ 
trix  (a  vector  containing  the  diagonal) ,  the  D  matrix,  if  diago¬ 
nal  and  the  matrix  option  is  chosen,  and  the  initial  conditions 
u°  and  uP  ,  if  present,  are  read  from  mass  storage. 

The  subroutine  CENDIF,  called  by  STINTC,  contains  the  in¬ 
tegrator  formulas  for  the  fixed  time  increment  and  increasing  or 
decreasing  time  increments;  see  13,4)  for  the  details  of  the 
various  formulas.  Both  un  and  u1*-5*  are  computed  based  on  the 
acceleration  at  tn_1  .  The  acceleration  is  computed  in  subrou¬ 
tine  ACTLN  for  the  cases:  1)  no  damping,  2)  diagonal  damping, 
and  3)  nondiagonal  damping;  again,  the  details  can  be  found  in 
references.  13,4).  Before  the  new  displacement'  un  and  the  velo- 

•  n-*5 

city  are  accepted,  the  error  for  this  step  is  computed  in 
subroutine  ERRORE  which  is  discussed  below.  Based  on  the  error 
computation  the  time  increment  for  the  next  step  is  computed  in 
subroutine  STEPSZ;  see  (4)  for  the  details.  Once  the  step  is 
accepted  the  computed  *  displacement  un  and  velocity 
in  (extrapolated  to  the  same  time  as  the  displacement)  and 
tn  are  transferred  through  the  user  written  subroutine  OUTPUT 
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for  display  by  printing,  plotting,  etc. 


In  subroutine  ERRORE  the  "maximum  perturbed  apparent”  fre- 
v  quency  C3,4J  and  the  "dominant  apparent  frequency*  (a  new  fre¬ 
quency  measure,  discussed  below)  error  measures  are  computed. 
From  14)  the  "maximum  perturbed  apparent  frequency"  error  meas¬ 
ure  is 
m 


e  =  max  (el1 
m  l 


:  MAXDEG 


) 


(2) 


where. 


€i  ”  <hn/4)  (ai/bi) 


aj  ”  I  a"  -  “T1  I 


. n  max 

bi=  u-j' 


-mb 


,  i  n  n-1  i . 
(I  u.  -u.  I) 


and  MAXDEG  is  the  size  (dimension)  of  the  solution  vector  in 
eqn.  (1),  hn  is  the  n-th  time  increment,  and  mb  is  the 
bandwidth  (average)  of  the  i-th  degree  of  freedom.  The  "domi¬ 
nant  apparent  frequency"  error  measure  £  a  is 


V- 


(A  un) 


(A  un) 


T 

T 


M  A  U 
**  — n 

M  A  un 


(3) 


where 
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A  un  =  un  -  u”"1 

and  superscript  T  indicates  transposition.  The  argument  of  the 
square  root  in  eqn.  (3)  is  recognized  as  a  form  of  Rayliegh's 
quotient.  Hence  it  gives  an  estimate  of  the  frequency  of  the 
global  (dominant)  change  in  the  solution  vector  un  at  tn  .  This 
measure  has  been  found  to  accurately  estimate  the  dominant  fre¬ 
quency  of  the  response.  The  user  controls  the  stability  and  ac¬ 
curacy  by  specifying  the  number  of  samples/cycle  for  both  the 
maximum  and  dominant  apparent  frequencies;  see  [3,4]  for  the 
relationship  between  samples/cycle  and  the  error  measures  given 
by  eqns.  (2)  and  (3).  Suggested  values  for  the  samples/cycle 
are  given  in  the  listing  for  the  user  written  subroutine  DRIVER 
and  the  integrator  performance  for  various  samples/cycles  is 
presented  with  the  sample  problem  in  Section  4. 

The  two  other  subroutines  ROTATE  and  VIPDA  are  utility  su¬ 
broutines  that  update  vector  address  pointers  and  compute  the 
vector  inner  product,  respectively. 

The  subroutines  listed  in  the  Appendices  are  the  DEC 
VAX-11/780  versions.  The  only  changes  required  to  use  these  su¬ 
broutines  on  a  UNIVAC  or  CDC  computer  are:  1)  the  INCLUDE  com¬ 
mand  for  inserting  the  labelled  common  CENTDF  into  the  subrou¬ 
tines,  and  2)  the  inline  comment  symbol  !  in  subroutine  STEPSZ 
and  the  user  written  subroutines  (Appendix  B) . 
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Limitations 

A  fixed  or  variable  time  increment  may  be  selected.  For 
|  the  fixed  time  increment  no  check  is  made  on  the  accuracy  or 
stability.  For  the  variable  time  increment  the  user  specifies 
the  minimum  and  maximum  time  increment  and  samples/cycle  for  the 
maximum  and  dominant  response  frequencies.  It  is  recommended 
that  the  minimum  time  increment  be  a  conservative  estimate  (very 
small) ,  because  the  automatic  time  increment  strategy  will  work 
more  effectively  with  little  loss  in  efficiency.  If  more  runs 
are  made  on  the  same  or  similar  problems,  the  minimum  time  in¬ 
crement  can  be  increased  if  the  performance  of  the  previous  run 
indicates  the  initial  choice  was  too  conservative. 

This  algorithm  will  constrain  a  degree  of  freedom  to  a  zero 
value  if  the  corresponding  mass  matrix  entry  is  zero.  No  other 
constraints  are  allowed. 

The  diagonal  damping  option  has  received  only  minimal  usage 
so  users  should  study  results  from  this  option  before  accepting 
them.  The  no  damping  and  nondiagonal  damping  options  have  been 
exercised  for  many  problems  and  the  authors  have  confidence  in 
their  performance. 

If  the  algorithm  is  used  in  the  variable  time  increment 
mode  for  pure  wave  propagation  response  (i.e.  integration  of 
the  linear  wave  equation) ,  it  will  generally  not  perform  at  its 
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best.  in  some  cases  the  results  are  terrible;  see  141  for  an 
example.  Note  that,  the  performance  for  this  problem  with  a 
properly  chosen  fixed  time  increment  is  excellent. 


3.  User  Written  Subroutines 


The  user  written  subroutines  are:  DRIVER,  FORCE,  DFORCE, 
S FORCE,  and  OUTPUT;  they  are  shown  schemetically  in  Figures  1 
and  2.  First,  DRIVER  transmits  the  problem  parameters  to  the 
time  integrator  through  the  arrays  INTGR,  LOGIC,  and  REALN,  plus 
the  starting  address  in  core  that  is  to*  be  used  by  the  time  in¬ 
tegrator.  Second,  the  applied  load  subroutine  FORCE,  provides 
the  load  vector  data  to  the  time  integrator.  Third,  subroutine 
DFORCE  provides  the  damping  force,  and  fourth,  subroutine  SFORCE 
provides  the  stiffness  force.  Finally,  subroutine  OUTPUT  pro¬ 
vides  display  of  the  response.  Note  that,  all  data  is  transmit¬ 
ted  as  vectors,  inasmuch  as  diagonal  matrices  may  be  considered 
vectors. 

The  listing  of  the  user  written  subroutines  in  Appendix  B 
includes  the  purpose  and  usage  requirements  of  each  subroutine. 
In  this  case  the  DRIVER  subroutine  also  includes  the  mass  matrix 
and  initial  velocity  data  for  the  example  problem:  normally 
this  data  is  generated  in  the  host  structural  analyzer.  Here 
the  structural  analyzer  is  embedded  in  the  user  written  inter- 
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face  subroutines  and  the  DRIVER  subroutine  is  the  main  program. 
In  a  typical  application  the  user  written  subroutine  DRIVER  is 
called  by  the  host  structural  analyzer  and  the  user  written  su- 
broutines  FORCE,  D FORCE ,  SFORCE,  and  OUTPUT  call  thie  host  struc¬ 
tural  analyzer  to  furnish  the  required  data. 


4.  Sample  Problem 


The  sample  problem  is  the  two  degree  of  freedom  nonlinear 
model  shown  in  Figure  3.  For  this  model  the  kinetic  energy,  T, 
potential  energy,  V,  and  the  dissipation  function,  F,  are  given 
by 


T  =  h  +  *5m2  Uj 

2 

V  ■  log  (cosh  u^)  +  ^  kc  -  u2)  +  k2  cosh  u2  (4) 


F  =  \  d  (Uj  -  u2) 

From  Lagrange's  equation  (5], 


the  mass  matrix,  damping  forces,  and  stiffness  forces  are  found 
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to  be 


dl  =  d^i 


*2> 


d2  =  d(u2  -  u1) 


(6) 


f,  =  k,  tanh  u,  +  k  (u,  -  u0) 

11  1  c  l  2 

f2  =  k2  cosh  u2  +  kc  (u2  -  u^) 

where  012  =  1112=  1.0,  d  ■  5.0,  “  k2  ~  1°°0*°  and  kc  =  100.0.  The 

initial  conditions  and  the  forcing  function  are  shown  in  Figure 
3. 


Physically,  the  system  comprises  a  softening  element  with 
an  initial  velocity  and  a  hardening  element  subjected  to  a  time 
delayed  rectangular  force  with  the  two  nonlinear  elements  moder¬ 
ately  coupled  by  a  linear  spring  and  dashpot.  This  model  was 
chosen  because  the  response  illustrates  how  well  the  time  incre¬ 
ment  selection  strategy  of  the  time  integrator  package  works. 

The  printed  output  for  the  sample  problem  is  listed  in  Ap¬ 
pendix  C.  The  first  three  lines,  the  time  increment  data  and 
the  synopsis  data  are  printed  from  the  STINT/CD  package.  The 
response  data  are  printed  in  subroutine  OUTPUT.  The  example 
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problem  was  also  run  for  20,  100,  and  500  samples/cycle  on  the 
dominant  frequency  (REALN( 5) )  in  addition  to  the  50 
saraples/cycles  run  shown  in  Appendix  C.  All  other  problem  par- 
?;  ameters  are  held  constant.  In  Figure  4  the  time  increment 
versus  the  response  time  is  shown.  In  general  the  time  incre¬ 
ment  is  increased  at  the  beginning  followed  by  decreases  after 
t  *  0.5  sec.,  when  the  rectangular  load  is  applied  to  the  har¬ 
dening  element.  This  behavior  is  very  reasonable.  A  large  time 
increment  is  possible  for  t  <  0.5  sec.  because  the  softening 
element  dominates  the  response,  but  later  the  hardening  element, 
with  a  higher  frequency,  dominates  the  response.  Also,  note 
that  as  the  sampling  rate  is  increased  (more  accuracy)  the  time 
increment  consistently  decreases. 

To  illustrate  that  convergence  is  obtained  the  example 
problem  was  also  run  for  a  fixed  time  increment  of  0.005  sec. 
This  time  increment  is  much  smaller  than  any  time  increment  se¬ 
lected  by  the  time  integrator  algorithm,  so  it  should  be  suffi¬ 
ciently  small  to  provide  a  converged  solution.  The  displacement 
and  velocity  history  for  the  four  variable  time  increment  runs 
and  the  fixed  time  increment  run  are  shown  in  Figures  5-8. 
During  the  first  half  of  the  response  the  results  are  nearly 
identical,  but  afterward  some  slight  differences  are  seen.  Note 
that  the  response  for  the  500  samples/cycle  and  the  fixed  time 
increment  are  identical;  look  near  the  very  end  of  the  response 
histories  for  the  clearest  picture. 
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In  a  fixed  time  increment  mode  this  problem  will  remain 
stable  for  a  time  increment  up  to  approximately  0.03  sec.  At 
this  time  increment  though  the  accuracy  is  minimal.  For  a  non¬ 
stiff  (closely  spaced  eigenvalues)  problem  with  a  small  (1-3) 
number  of  degrees  of  freedom  accuracy  considerations  usually 
dominate  over  stability  considerations.  For  a  stiff  (widely 
separated  eigenvalues)  problem  with  more  degrees  of  freedom, 
such  as  the  cantilever  beam  [4J,  stability  dominates  and  the  ac¬ 
curacy  is  well  within  what  the  user  requests. 

For  the  variable  time  increment  mode  the  behavior  of  the 
average  time  increment  versus  the  requested  samples/cycle  of  the 
dominant  apparent  frequency,  shown  in  the  Table  below,  sheds 
more  light  on  stability  versus  accuracy. 


Dominant  Frequency 
Samples/Cycle 

20 

50 

100 

500 


Average  Time 
Increment  (sec.) 

0.0062112 

0.0056818 

0.0045455 

0.0010917 


Only  the  time  increment  change  from  100  to  500  samples/cycle 
shows  a  decrease  in  proportion  to  the  samples/cycle  ratio.  At 
the  other  sampling  rates  the  time  increment  is  being  controlled 
by  a  mixture  of  stability  and  accuracy.  Hence  the  ratios 
between  the  sampling  rate  and  the  average  time  increment  do  not 
correlate  unless  accuracy  clearly  controls  the  time  increment. 
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5.  Suggested  Usage 


fe;  i, 

v  The  stand-alone  in-core  time  integrator  package  STINT/CD 

presented  here  is  easily  adapted  to  problems  with  several  hun¬ 
dred  degrees  of  freedom.  However,  before  attempting  a  problem 
this  large,  it  is  suggested  that  the  user  experiments  with  a 
small  problem  with  characteristics  similar  to  a  larger  problem 
that  one  is  interested  in  solving.  The  user  may  also  want  to 
change  some  of  the  algorithm  parameters  in  the  subroutine  STEPSZ 
to  fit  a  particular  class  of  problems  more  efficiently. 
References  [3,4]  should  be  studied  before  adjusting  the  time  in¬ 
crement  selection  parameters. 

In  the  authors'  work  environment  (solving  a  variety  of  in¬ 
terdisciplinary  problems)  the  stand-alone  feature  is  most  de- 
sireable,  since  it  allows  the  quick  assembly  of  software  ele¬ 
ments  to  solve  the  current  problem.  However,  in  a  strictly  pro¬ 
duction  environment,  where  the  mechanics  of  the  problem  seldom 
change,  some  efficiency  may  be  gained  by  embedding  the  time  in¬ 
tegrator  into  the  structural  analyzer;  see  (61. 

The  explicit  central  difference  integrator  is  most  suited 
to  short  duration  transient  loads  or  problems  in  which  the 
stiffness  may  suddenly  change.  The  automatic  variable  time  in¬ 
crement  feature  is  especially  suited  to  these  problems  in  that 
it  selects  the  largest  possible  time  increment  consistent  with 
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maintaining  stability  and  accuracy  based  on  the  current  system 
behavior.  Therefore  a  small  time  increment  does  not  have  to  be 
used  during  the  entire  computations  to  achieve  accuracy  during 
one  small  time  interval,  since  the  small  time  increment  is  auto¬ 
matically  selected  only  when  it  is  needed. 
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Appendix  A 

STINT/CD  Computer  Code  Listing 


SUBROUTINE  STINTC(LOGIC, INTGR, REALN,C) 

***  PURPOSE  PROVIDES  THE  INTERFACE  BETWEEN  THE  USER  SUPPLIED 

'DRIVER'  ROUTINE  AND  THE  CENTRAL  DIFFERENCE  TIME 
INTEGRATOR  CENDIF/VECTOR  (CALL  CENDIF) 

***  INPUT  LOGIC  ARRAY  OF  LOGICAL  CONTROL  DATA 

INTGR  ARRAY  OF  INTEGER  CONTROL  DATA 
REALN  ARRAY  OF  REAL  CONTROL  DATA 

C  STARTING  ADDRESS  OF  COMMON  (CORE)  TO  BE  USED  BY 
THE  INTEGRATOR 


EXTERNAL  GASPER 
INCLUDE  ' PRODCD. PCR ' 


LOGICAL  LOGIC (1) 

LOGICAL  IGASP, IDISP, IVEL 
INTEGER  INTGR(l) 

REAL  C(l) 

REAL  REALN ( 1 ) 

C 

IGASP  «  LOGIC (1) 

I FORCE  «  LOGIC (2) 

IDISP  «  LOGIC (3) 

IVEL  -  LOGIC (4) 

FIXSTP  «  LOGIC (5) 

IDAMP  ■  LOGIC (6) 

IDMPDG  ■  LOGIC (7) 

C 

MAXDEG  ■  INTGR(l) 

IUNIT  ■  INTGR(2) 

KBAND  ■  INTGR<3) 

KBAND  «  KBAND  -  1 

C 

TIME  ■  REALN ( 1 ) 

TMAX  ■  REALN (2) 


non 
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DTMIN 

DTMAX 

EPSDFQ 

EPSHFQ 

ALFA 


-  REALN( 3) 
*  REALNC4) 
)  =  REALN(5) 
)  *  REALN(6) 
«  REALN ( 7 ) 
0.5* ALFA 


IF (FIXSTP)  DTMIN  *  DTMAX 
DT  *  DTMIN 

MAXSTP  *  3*TMAX/DTMIN 
NSTEP  =  1 

IF (EPSHFQ  .LT.  3.14)  EPSHFQ  *  4.0 
EPSHFQ  =  (3.14159265/EPSHFQ) **2 
IF (EPSDFQ  .LT.  3.14)  EPSDFQ  *  4.0 
EPSDFQ  *  (3 . 141 5926 5/EPSDFQ) **2 
UPNTR(l)  =  0 
UPNTR( 2)  -  1 
OPNTR ( 3 )  »  2 


UPNTR( 2)  ■ 
OPNTR ( 3 )  * 
MAXVEC  *  3 
IF(FIXSTP) 
IF(FIXSTP) 
IF(FIXSTP) 
NDOUBL  -  0 
NCOTS  *  0 
NFACTS  -  0 


MAXVEC  *  1 
OPNTR (2)  = 
OPNTR (3)  * 


***  STORAGE  ALLOCATION 

M  *  MAXDEG 

MM  *  MAXVEC *MAXDEG 

LSM  =  1 

LRSM  *  LSM  +  M 
IF (FIXSTP)  LRSM  *  LSM 
LSC  ■  LRSM  +  M 
LU  *  LSC 

IF(IDMPDG)  LO  ■  LSC  +  M 
LUD  *  LU  +  MM 
LUDD  »  LUD  +  MM 
LW  *  LUDD  +  MM 
LWORKA  »  LW 

IF( (IDMPDG)  .OR.  (IDAMP))  LWORKA  *  LW  +  M 
LWORKB  «  LWORKA  +  M 
LEND  -  LWORKB  +  M  -  1 

WRITE (6, 200)  LEND 

200  FORMAT (1H1 ,  2X,’WELCOME  TO  STINTC  THE  STAND-ALONE  CENTRAL  DIFFEREN 
ICE  TIME  INTEGRATOR' ,//,17X, 'STINTC  REQUIRES' , 110, '  WORDS  OF  STORA 
2GE' ,//,27X, '6  DECEMBER  1979  VERSION',//) 

DO  10  1*1 , LEND 
C(I>  *  0.0 
10  CONTINUE 
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IF(IGASP)  GO  TO  20 

C  ***  UNFORMATTED  READS 

READ (I UNIT)  (C(LSH-1+I> ,1-1, M) 

IF ( (IDMPDG)  .AND.  (IDAMP))  READ(IUNIT)  (C (LSC-l+I) , 1-1 ,M)  i 

IF(IDISP)  READ(IUNIT)  (C(LU“1+I) , I»I ,M) 

IF(IVEL)  READ (I UNIT)  (C (LUD-l+I) , 1=1 ,M)  k  j 

CLOSE ( UNIT- IUNIT) 

GO  TO  30 

***  DMGASP  READS  NOTE:  IF  NOT  AVAILABLE  REPLACE  BY  OTHER  EFFICIENT 
***  I/O  PACKAGE 

20  CALL  DM  RAST  <IUNIT,  C(LSM) ,  M) 

CALL  ETS  PRT ( 1 » 1 STINTCR1 * ) 

CALL  DM  HAST  (0,  GASPER,  0) 

IF( (IDMPDG)  .AND.  (IDAMP))  CALL  DM  RAST  (IUNIT,  C(LSC),  M) 

CALL  ETS  PRT ( 1 , 1 STINTCR2 ' ) 

CALL  DM  HAST  (0,  GASPER,  0) 

IF(IDISP)  CALL  DM  RAST  (IUNIT,  C(LU) ,  M) 

CALL  ETS  PRT ( 1 , 1 STINTCR3  * ) 

CALL  DM  HAST  (0,  GASPER,  0) 

IF(IVEL)  CALL  DM  RAST  (IUNIT,  C(LUD),  M) 

CALL  ETS  PRT(1, 'STINTCR4 ') 

CALL  DM  HAST  (0,  GASPER,  0) 

CALL  DM  FAST  (IUNIT,  0,  0) 

30  CONTINUE 


CALL  CENDIF(  C(LSM) ,  C(LRSM),  C(LSC),  C(LU),  C(LUD) ,  C(LUDD), 
1  C ( LW) ,  C ( LWORKA) ,  C ( LWORKB ) ) 

RETURN 

END 


SUBROUTINE  CENDIF ( SM,  RSM,  SC,  U,  UD,  UDD,  W,  WORKA,  WO.lIvB) 


SOLVES  SECOND  ORDER  STRUCTURAL  DYNAMICS  BY  EXPLICIT  CENTRAL 
DIFFERENCE  METHOD 


***  EQUATIONS  OF  MOTION 

M*UDD  +  D*UD  +  K*U  «  F 

WHERE  M  IS  A  DIAGONAL  MASS  MATRIX 
D  IS  THE  DAMPING  MATRIX 

K*U  IS  THE  FORCE  DUE  TO  STIFFNESS  (LINEAR  OR  NONLINEAR)! 
***  VECTOR  FORMULATION  (IE  D*UD  AND  K*U  ARE  AVAILABLE  AS  VECTORS) 

***  DEFINITIONS 

W  -  UD  +  0 ,5*ALFA* (HN+HNM1) * (DV  -  D*UD/M)  (NONDIAGONAL 
DAMPING) 
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W  USED  TO  STORE  DAMPING  FORCE  (DIAGONAL  DAMPING) 

HN  THE  TIME  STEP  FROM  N  TO  N+HALF 
HNH1  THE  TIME  STEP  FROM  N-HALF  TO  N 

***  STORAGE 

SM  MASS  MATRIX  (DIAGONAL  STORED  AS  A  VECTOR) 

RSM  RECIPROCAL  OF  MASS  MATRIX 

SC  DAMPING  MATRIX  (ONLY  REQUIRED  FOR  DIAGONAL  DAMPING,  STORED 
AS  A  VECTOR) 

U  (3  MOST  CURRENT  VALUES  FOR  VARIABLE  STEP) 

UD  (3  MOST  CURRENT  VALUES  FOR  VARIABLE  STEP) 

UDD  (3  MOST  CURRENT  VALUES  FOR  VARIABLE  STEP) 

W  SEE  DEFINITION  ABOVE 
SCRATCH  VECTORS  WORKA, WORKB 

INCLUDE  'PRODCD.PCR' 


INTEGER  NTIME(28) 

REAL  SM ( 1 ) 

REAL  RSM ( 1 ) 

REAL  SC (I) 

REAL  U(l) , UD ( 1 ) , UDD ( 1 ) 

REAL  W(l) 

REAL  WORKA (1) , WORKB (1) 

***  initial  VALUES 
NSOLV  «  0 
IRST  »  0 
I LAST  -  0 
IDTMX  »  0 
INCTRY  «  0 
IDEC  ■  0 
INCRS  »  .FALSE. 

DECRS  »  .FALSE. 

ERR  -  0.0 
RN  «  1. 

C  ***  INITIALIZE  TIME  STEP  OCCURENCE  ARRAY 
DO  10  1-1,28 
NTIME(I)  -  0 
10  CONTINUE 

C  ***  FORM  INVERSE  OF  DIAGONAL  MASS  MATRIX  SM 
C  ***  NOTE  IF  SM-0.0  DOF  IS  CONSTRAINED  TO  0.0 
DO  20  1-1 , MAXDEG 

IF(SM(I)  .EQ.  0.0)  GO  TO  20 

RSM ( I )  -  l./SM(I) 

20  CONTINUE 
C  ***  STEP  SIZES 
30  CONTINUE 

IF(IRST.EQ.O)  HN  -  5.*DTMIN 
IF(HN  ,.GT.  (0.5*DTMAX) )  HN  ■  0.5*DTMAX 
IF(FIXSTP)  HN  -  0.5*DTMAX 
HNM1  -  HN 
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***  STARTING  PROCEDURE 


***  INITIALIZE  WORK  SPACE  . 

DO  100  1*1 , MAXDEG  f 

WORK A (I)  *0.0 
WORKB(I)  *  0.0 

***  W  *  1  FOR  DIAGONAL  DAMPING  INITIAL  SET  UP  (IF  IDAMP  FALSE) 
W(I)  *  1.0 
100  CONTINUE 

***  COMPUTE  APPLIED  FORCE  IN  WORKA 

IF ( IFORCE)  CALL  FORCE (MAXDEG, TIME, WORKA) 

***  COMPUTE  TOTAL  SPRING  FORCE  (IN  WORKB) 

CALL  SFORCE (MAXDEG, U, WORKB) 

***  PRINT  OUT  INITIAL  CONDITIONS 
I PRT* 2 

CALL  OUTPUT ( I PRT, MAXDEG , TIME, U, UD) 

IF(IRST  .EQ.  0)  CALL  OUTPUT (I PRT, MAXDEG, TIME, U,UD) 

***  FOR  DIAGONAL  DAMPING  PUT  DAMPING  COEFFICIENTS  INTO  SC 

IF ( ( IDMPDG)  .AND.  (.NOT.  IDAMP))  CALL  DFORCE (MAXDEG, W, SC) 

***  COMPUTE  DAMPING  FORCE  (IN  W) 

IF( (IDAMP)  .AND.  (.NOT.  IDMPDG))  CALL  DFORCE (MAXDEG, UD,W) 

***  FORM  -F+D*DU+K*U  IN  WORKA 
DO  110  1*1, MAXDEG 

IF (.NOT. IDAMP)  W(I)  *  0.0 
IF ( IDMPDG)  W(I)  -  SC (I) *UD(I) 

WORKA(I)  *  -WORKA(I)  +WORKB(I)  +W(I) 

110  CONTINUE 

***  COMPUTE  ACCELERATION  AT  TIME* START,  VELOCITY  AT  TIME* ST ART+  H N , 
***  AND  DISPLACEMENT  AT  TIME*START+2HN 

CALL  ROTATE ( MAXVEC , UPNTR) 

11  *  UPNTR ( 1) *MAXDEG 

12  *  UPNTR ( 2 ) * MAXDEG 
DO  120  1*1, MAXDEG 

J*I1+I 

K*I2+I 

UDD(K)  *  -RSM ( I ) *WORK A ( I ) 

UD(J)  *  UD (K)  +  HN*UDD(K) 

U(J)  *  U (K)  +  2 . *HN*UD ( J) 

120  CONTINUE 

***  BEGIN  TIME  LOOP 

N STEPS  -  NSTEP-1 
200  NSTEPS  *  NSTEPS+1 


***  ARRAY  POINTERS 

11  -  UPNTR ( 1 ) * MAXDEG 
III  *  .11+1 

12  ■  UPNTR ( 2 ) * MAXDEG 

13  -  UPNTR ( 3) +MAXDEG 
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C  ***  UPDATE  TIME 

300  TIME  »  TIME  +  2.*HN 
NSOLV  «  NSOLV  +  1 
IF (NSTEPS  .EQ.  NSTEP) 


GO  TO  340 


***  UPDATE  VELOCITY  AND  DISPLACEMENT 

***  (VELOCITY  AT  N-HALF,  DISPLACEMENT  AT  N) 

DO  330  Is 1 r MAXDEG 
J«I1+I 
K*I2+I 
L*I3+I 

IF(RN  .LE.  1.0)  GO  TO  310 

D  =  ABS(UDD(K) ) 

IF (D  .EQ.  0.0)  GO  TO  310 

DEM  *  ABS(UDD(K) -UDD(L) ) /D 
IF (DEM  .GT.  l.E-02)  GO  TO  310 

***  ALMOST  CONSTANT  ACCELERATION 

UD ( J)  »  UD(K)-0.25*(HN-HNMl)*((3.-RN)*UDD(K)+(l.+RN)*UDD(L) ) 
UD ( J)  =  UD(J)  +  2.0*HN*UDD(K) 

GO  TO  320 

310  UD ( J)  *  UD (K)  +  (HN+HNM1) *UDD(K) 

320  U(J>  *  U(K)  +  2 . 0*HN*UD ( J) 

330  CONTINUE 
340  CONTINUE 

***  COMPUTE  ACCELERATION 

CALL  ACLTNlRSM, SC, U, UD, UDD , W, WORKA, WORKB) 

***  ERROR  COMPUTATION 

IF (FIXSTP)  GO  TO  430 

CALL  ERRORE(SM, U, UDD, ETRNM, ETRDF, WORKA, WORKB) 

***  STEP  SIZE  CONTROL  SECTION 

CALL  STEPS Z ( ETRNM, ETRDF, Uf  UD, UDD,KGO) 

GOTO  (400,30,500,300) ,KGO 


***  OUTPUT  SECTION 


400  CONTINUE 
IDEC  »  0 

C  ***  TIME  STEP  DISTRIBUTION  DETERMINATION 
DT  -  2 .0*HNM1 
TLIM  »  0.0 
TFAC  -  10.0 
TLIMA  -1.0 
DO  410  INT-1 ,28 

IF ( INT  .GT.  10)  TFAC  *  100. 

IF(INT  .GT.  19)  TFAC  -  1000. 

TLIM  -  TLIM  +  TFAC 

IF ( (DT  ,GE.  TLIMA*DTMIN)  .AND.  (DT  .LT.  TLItl*DTMIN)  ) 

GO  TO  420 
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TLIMA  -  TLIM 
410  CONTINUE 

420  NTIME(INT)  *  NTIME(INT)  +  1 
430  CONTINUE 
IPRT  «  0 

£  *★*  SET  TIME  STEP  TO  END  AT  TMAX 
'  IF (TIME  .GE.  TMAX)  I LAST* 1 

IF( (TIME+2 .0*HN)  .GT.  TMAX)  HN  =  0 . 5* (TMAX-TIME) 
C  ***  PRINTING  OUTPUT 

C  ***  COMPUTE  VELOCITY  AT  N-STEP  FOR  PRINTOUT  PURPOSES 
DO  440  I®1 , MAXDEG 
J*I1+I 

WORKB(I)  =  UD ( J)  +  HN*UDD ( J) 

440  CONTINUE 

IF ( ILAST  .EQ.  1)  IPRT® 10 

CALL  OUTPUT ( I PRT , MAXDEG , TIME ,  U  ( 1 1 1 ) ,WORKB) 

IF (ILAST  .EQ.  1)  GO  TO  500 

C  ***  SET  POINTERS  FOR  NEXT  TIME  STEP 
CALL  ROTATE (MAXVEC, UPNTR) 


***  END  TIME  LOOP 

IF (NSTEPS  .LT.  MAXSTP)  GO  TO  200 

***  SUMMARY  PRINT  OUT  OF  TIME  INTEGRATOR  PERFORMANCE 
500  IF(FIXSTP)  GO  TO  520 

TAVE  »  TIME/ (NSTEPS-1) 

WRITE(6,1000) 

1000  FORMAT (1 HI) 

WRITE (6,1010)  TAVE , NDOUBL , NCUTS , NFACTS , NSTEPS , NSOLV 


1010  FORMATC  AVERAGE  TIME  STEP  ',E12.5/, 

12X,  'NUMBER  OF  STEP  INCREASES  'rI5/, 

22X,  'NUMBER  OF  STEP  DECREASES  ',15/, 

32X,  'NUMBER  OF  FACTORIZATIONS  ',15/, 

42X,  'NUMBER  OF  TIME  STEPS  ',15/, 

52X,  'NUMBER  OF  SOLUTIONS  ’,15  ) 

WRITE (6, 1020) 


1020  FORMAT (////, 1 OX, 'DT  OCCURRENCES  IN  THE  RANGES  INDICATED',/) 

TLIM  =  0.0 
TFAC  *  10.0 
TLIMA  «  1.0 
DO  510  INTIM*1 ,28 

IF ( INTIM  .GT.  10)  TFAC  =  100. 

IF ( INTI M  .GT.  19)  TFAC  *  1000. 

TLIM  «  TLIM  +  TFAC 
NTI  ®  NTIME( INTIM) 

IF (NTI  .NE.  0)  WRITE (6, 1030)  TLIMA, TLIM, NTI 
1030  FORMATC  FROM' , F8 . 1 ,2X ,' TO' , F8 . 1 ,2X, ' TIMES  DTMIN,  DT  OCCURREMC 
1ES  WAS' ,15) 

TLIMA  ■  TLIM 
510  CONTINUE 
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520  RETURN 
END 


SUBROUTINE  ACLTN ( RSM , SC , U , UD , UDD , W , WORK A, WORKB) 

***  PURPOSE  TO  COMPUTE  THE  ACCELERATION  AT  THE  N-TH  STEP 

***  INPUT  HN  HALF  OF  CURRENT  TIME  STEP  (COMMON  CENTDF) 

HNM1  HALF  OF  PAST  TIME  STEP  (COMMON  CENTDF) 
TIME  CURRENT  TIME  (COMMON  CENTDF) 

RSM  RECIPROCAL  MASS  MATRIX  ARRAY 
SC  DIAGONAL  DAMPING  MATRIX  ARRAY 
U  DISPLACEMENT  ARRAY 
UD  VELOCITY  ARRAY 

***  SCRATCH  W  ARRAY  FOR  DAMPING 

WORK A, WORKB 

***  OUTPUT  UDD  ACCELERATION  ARRAY 

INCLUDE  'PRODCD.PCR' 


REAL  RSM(l) ,SC(1) 

REAL  U(l) ,  UD ( 1 ) , UDD ( 1 ) 

REAL  W(l) 

REAL  WORKA <  1 )  , WORKB ( 1 ) 

C 

III  -  11+1 
C 

C  ***  INITIALIZE  WORKING  ARRAYS  (PREPARING  TO  COMPUTE  ACCELERATION 
C  ***  AT  N-TH  STEP) 

DO  10  1=1 , MAXDEG 
WORKA(I)  »  0.0 
WORKB ( I )  =  0.0 
10  CONTINUE 

C  ***  COMPUTE  APPLIED  FORCE  IN  WORKA 

IF ( I FORCE)  CALL  FORCE (MAXDEG, TIME, WORKA) 

C  ***  COMPUTE  TOTAL  SPRING  FORCE  IN  WORKB 
CALL  SFORCE (MAXDEG , U ( I 11 ) , WORKB) 

C  ***  FORM  -F+K*U  ( -WORK A+WORKB )  AND  FORM  DV  (IN  WORKA) 

DO  20  1=1, MAXDEG 

WORKA ( I)  =  -RSM ( I ) * (-WORKA ( I ) +WORKB ( I ) ) 

20  CONTINUE 

IF ( IDMPDG)  GO  TO  120 

IF( .NOT.  IDAMP)  GO  TO  100 

C  ***  COMPUTE  DAMPING  FORCE  IN  WORKB 

CALL  DFORCE (MAXDEG , UD (III) , WORKB) 

C  ***  FORM  W 

DO  30  1=1, MAXDEG 
J-Il+I 

W(I)  ■  RSM ( I ) *WORKB ( I ) 

W(I>  ■  UD ( J)  +  BTA* (HN+HNMl ) * (WORKA(I) -W( I ) ) 
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30  CONTINUE 

C  ***  COMPUTE  DAMPING  AS  UPDATED 
CALL  D FORCE ( MAXDEG , W , WORKB) 

GO  TO  150 

***  COMPUTE  ACCELERATION  AT  N 

***  NO  DAMPING 
100  DO  110  1=1, MAXDEG 
J-Il+I 

UDD(J)  =  WORK A ( I ) 

110  CONTINUE 

GO  TO  200 

***  DIAGONAL  DAMPING 
120  DO  130  1=1, MAXDEG 
J=I1+I 

FI  =  1.  +  HN*RSM(I) *SC(I) 

F2  =  RSM ( I ) *SC ( I ) 

WORKB ( I )  =  WORKA(I)-  F2*UD(J) 

WORKB ( I)  =  UD ( J) +  HN*WORKB(I)/Fl 
130  CONTINUE 

***  DETERMINE  ACCELERATION 

CALL  DFORCE (MAXDEG, WORKB, W> 

DO  140  1=1, MAXDEG 
J=I1+I 

UDD(J)  =  WORKA(I)  -  RSM ( I ) *W ( I ) 

140  CONTINUE 

GO  TO  200 

***  NONDIAGONAL  DAMPING 
150  DO  160  1=1, MAXDEG 
J»I1+I 

UDD(J)  =  WORKA(I)  -  RSM ( I ) *WORKB ( I ) 

160  CONTINUE 

200  RETURN 
END 


SUBROUTINE  ERRORE (SM, U, UDD, ETRNM, ETRDF, WORK A, WORKB) 

***  PURPOSE  TO  COMPUTE  THE  ERROR  BASED  ON  THE  HIGHEST 

APPARENT  FREQUENCY  CONCEPT  AND  THE 
DOMINANT  FREQUENCY 

***  INPUT  HN  HALF  OF  CURRENT  TIME  STEP  (COMMON  CENTDF) 

SM  MASS  MATRIX 
U  DISPLACEMENT  ARRAY 
UDD  ACCELEREATION  ARRAY 
COMMON  /CENTDF/ 

***  OUTPUT  ETRNM  THE  ERROR  (MAX  LOCAL  FREQUENCY) 

ERTDF  THE  ERROR  (DOMINANT  FREQUENCY) 

***  SCRATCH 


WORK A, WORKB 
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INCLUDE  * PRODCD. PCR ' 

C 

REAL  SM ( 1 ) ,U(1) , UDD(l) ,WORKA(l) ,WORKB(l) 

APPROPRIATE  VALUES  OF  EPSMAC  FOR  VARIOUS  COMPUTERS  ARE 
CDC  6000/7000  SERIES  7.11E-15  (SINGLE  PRECISION) 

IBM  360/370  SERIES  9.54E-07  (REALM  PRECISION) 

IBM  360/370  SERIES  2.22E-16  (REAL* 8  PRECISION) 

UNIVAC  1108/1110,  IBM  7094  1.49E-00  (SINGLE  PRECISION) 


C 


C 


C 


DATA  EPSMAC  /9.54E-07/ 


HT24  *  HN*HN 
ETRNM  *  0.0 


***  FIND  MAXIMUM  LOCAL  FREQUENCY 
DO  30  1=1 , MAXDEG 
J=I1+I 
K=I2+I 


1 


L=I3+I 

WORK A ( I )  =  U(J)  -  U(K) 
WORKB(I)  =  SM (I ) * (UDD ( J) 
DEM  =  ABS(U( J) -U(K) ) 

IF (DEM  .EQ.  0.0) 

IF ( (U ( J)  .NE.  0.0)  .AND. 


-  UDD(K) ) 

GO  TO  30 
(DEM/ABS(U( J) ) 
GO  TO  30 
GO  TO  20 


I F ( KBAND  .EQ.  0) 

KBEG  =  I-KBAND 

IF (KBEG  .LE.  0)  KBEG  =  1 

KEND  =  I+KBAND 

IF (KEND  .GT.  MAXDEG)  KEND  =  MAXDEG 
DO  10  KBE=KBEG, KEND 
JJ  =  Il+KBE 
KK  =  I2+KBE 
D  =  ABS(U( JJ) -U(KK) ) 

DEM  =  AMAX1 (DEM,D) 


10  CONTINUE 

20  CONTINUE 

ERR  =  HT24* (UDD( J) -UDD(K) ) /DEM 
IF (ABS(ERR) .GT.ABS (ETRNM) )  ETRNM  =  ERR 
30  CONTINUE 

***  FIND  DOMINANT  FREQUENCY 

BOT2  »  VI PD A (WORK A, WORK A, MAXDEG) 

TOP  =  VI PDA (WORK A, WORKB, MAXDEG) 

DO  40  1=1, MAXDEG 

WORKB(I)  =  SM(I) *WORKA(I) 


.LE.  EPSMAC)) 


40  CONTINUE 

BOT  =  VI PDA (WORK A, WORKB, MAXDEG) 

I F ( BOT2  .EQ.  0.0)  BOT  =  0.0 
I F ( BOT2  .EQ.  0.0)  BOT2  «  1.0 

IF( (ABS(BOT) /BOT2)  .LT.  ( FLOAT(MAXDEG) *EPSMAC) )  BOT  =  0. 
ALAMS  =0.0 

IF (BOT  .NE.  0.0)  ALAMS  =  TOP/BOT 
ETRDF  »  HT24*SQRT(ABS( ALAMS) ) 


C 
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30  CONTINUE 

C  ***  COMPUTE  DAMPING  AS  UPDATED 
CALL  D FORCE (MAXDEG, W, WORKB) 

GO  TO  150 

***  COMPUTE  ACCELERATION  AT  N 

***  NO  DAMPING 
ICO  DO  110  1=1, MAXDEG 
J=I1+I 

UDD(J)  =  WORK A ( I ) 

110  CONTINUE 

GO  TO  200 

***  DIAGONAL  DAMPING 
120  DO  130  1=1 , MAXDEG 
J=I1+I 

FI  =  1.  +  HN*RSM( I ) *SC (I ) 

F2  =  RSM ( I ) *SC ( I ) 

WORKB ( I )  =  WORKA(I)-  F2*UD(J) 

WORKB(I)  =  UD(J)+  HN*WORKB ( I ) /FI 
130  CONTINUE 

***  DETERMINE  ACCELERATION 

CALL  D FORCE ( MAXD EG, WORKB,W) 

DO  140  1=1 , MAXDEG 
J=I1+I 

UDD(J)  =  VZORKA(I)  -  RSM(I)  *W(I) 

140  CONTINUE 

GO  TO  200 

***  NONDIAGONAL  DAMPING 
150  DO  160  1=1, MAXDEG 
J=I1+I 

UDD(J)  =  WORKA(I)  -  RSM ( I ) *WORKB ( I ) 

160  CONTINUE 

200  RETURN 
END 


SUBROUTINE  ERRORE ( SM, U,UDD, ETRNM, ETRDF, WORKA, WORKB) 

***  PURPOSE  TO  COMPUTE  THE  ERROR  BASED  ON  THE  HIGHEST 

APPARENT  FREQUENCY  CONCEPT  AND  THE 
DOMINANT  FREQUENCY 

***  INPUT  HN  HALF  OF  CURRENT  TIME  STEP  (COMMON  CENTDF) 

SM  MASS  MATRIX 
U  DISPLACEMENT  ARRAY 
UDD  ACCELEREATION  ARRAY 
COMMON  /CENTDF/ 

***  OUTPUT  ETRNM  THE  ERROR  (MAX  LOCAL  FREQUENCY) 

ERTDF  THE  ERROR  (DOMINANT  FREQUENCY) 

***  SCRATCH 


WORKA, WORKB 
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RETURN 

END 


SUBROUTINE  STEPSZ (ETRNM, ETRDF, U, UD, UDD, KGO)  ' 

***  PURPOSE  STEPSIZE  SELECTION 

***  INPUT  HNM1  HALF  OF  PAST  TIME  STEP  (COMMON  CENTDF) 

TIME  CURRENT  TIME  (COMMON  CENTDF) 

ETRNM  CURRENT  ERROR  (HIGH  FREQUENCY) 

ETRDF  CURRENT  ERROR  (DOMINANT  FREQ) 

U  DISPLACEMENT  ARRAY 
UD  VELOCITY  ARRAY 
UDD  ACCELERATION  ARRAY 
COMMON  /CENTDF/ 

***  OUTPUT  HN  HALF  OF  NEW  TIME  STEP  (COMMON  CENTDF) 

RN  HN/HNMl  (COMMON  CENTDF) 

KGO  COMPUTED  GOTO  CONTROL 

INCLUDE  ' PRODCD. PCR' 

C 

REAL  U(l) , UD ( 1 ) , UDD ( 1 ) 

C 

TERM  *  ABS ( ETRNM) /EPSHFQ 
TERM2  =  ABS ( ETRDF) /EPSDFQ 
TERM  ■  AMAX1 (TERM, TERM2 ) 

C  ***  POSSIBLE  STEP  SIZE  INCREASE  CHECK 

IF (TERM  .LT.  0.1)  GO  TO  10 

C  ***  STEp  SIZE  DECREASE  CHECK 

IF (TERM  .GT.  1.0)  GO  TO  20 

HNM1  *  HN 

RN  =  1.0 

INCTRY  «  0 

IDEC  =  0 

INCRS  -  .FALSE. 

DECRS  *  .FALSE. 

GO  TO  50 

C  ***  STEP  SIZE  INCREASE 
10  HNM1  -  HN 

INCRS  »  .FALSE. 

INCTRY  *  INCTRY  +  1 
IDEC  **  0 

IF (NSTEPS  .LE.  6)  INCTRY  -  0 

C  ***  MUST  TRY  TO  INCREASE  FOR  5  CONSECUTIVE  STEPS  BEFORE  INCREASING 
IF ( INCTRY. LE. 4)  GO  TO  50 

INCTRY  =  0 
NDOUBL  *  NDOUBL  +  1 
INCRS  -  .TRUE. 

DECRS  =  .FALSE. 

IF (TERM  . EQ.  0.0)  TERM  ■  0.001 
RATIO  -  (1. 0/(4. *TERM) ) **0.25 
RATIO  -  AMIN1 (RATIO ,1.5) 
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C 


C 


C 


C 


HN  »  RATIO*HN 
RN  *  HN/HNM1 
DTN  *  2  .  *HN 

IF (DTN  .GT.  DTMAX)  HN  *  0.5*DTMAX 

IF (DTN  .GT.  DTMAX)  IDTMX  =  IDTMX+1 

IFdDTMX  .GT.  5)  GO  TO  50 

DTN  *  2 . *HN 

WRITE (6 ,200)  DTN, TIME 

200  FORMAT ( *  $$$  STEP  SIZE  INCREASED  TO  ',E10.4,' 

GO  TO  50 

***  STEP  SIZE  DECREASE 
20  TIME  *  TIME  -  2.0*HN 
INCTRY  =  0 
INCRS  *  .FALSE. 

DECRS  =  .TRUE. 

IDTMX  =  0 
IDEC  =  IDEC  +  I 

***  IF  DECREASE  MORE  THAN  4  TIMES  AT  ONE  TIME  POINT, 
IF (IDEC  .GT.  4)  GO  TO  30 

RATIO  *  (1.0/(5.*TERM) ) **0.2 
IF (RATIO  .LT.  0.66666667)  RATIO  =  0.66666667 
IF (RATIO  .GT.  0.9)  RATIO  =  0.9 
HN  =  RATIO*HN 
RN  *  HN/HNM1 
DTN  «  2 . *HN 
NCUTS  =  NCUTS  +  1 
WRITE (6,210)  DTN, TIME 

210  FORMAT ( '  $$$  STEP  SIZE  DECREASED  TO  ’,E10.4,’ 

IF (DTN  .LT.  DTMIN)  GO  TO  70 

IF (NSTEPS  .NE.  NSTEP)  GO  TO  80 

***  FIRST  TIME  STEP  SO  WE  ARE  GOING  TO  RESTART 
UPNTR(l)  =  0 
UPNTR ( 2 )  =  1 
UPNTR( 3)  =  2 
HNMl  =  HN 
IRST  =  1 
RN  =1.0 
WRITE (6,220) 

220  FORMAT (2X, *$$$', 2X, 'RESTART') 

GO  TO  €0 

30  CONTINUE 

***  LOAD  U  AND  UD  FOR  A  RESTART  AT  TIME  *  TIME 
DO  40  1=1 , MAXDEG 
K=I2+I 
U(I)  =  U(K) 

UD ( I)  =  UD (K)  +  HNM1*UDD(K) 

40  CONTINUE 

***  INITIALIZE  POINTERS 
UPNTR(l)  *  0 
UPNTR ( 2 )  =  1 
UPNTR ( 3)  =  2 
WRITE(6,230)  TIME 

230  FORM AT ( 2X,' $$$' /2X, ’RESTART  AT  TIME  »  ’,E12.5) 
IRST  ■  1 
IDEC  »  0 


AT  ' , El 0 . 4 ) 


RESTART 


AT  ' , E10 .4) 
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ERR  *  0.0 
NSTEP  *  NSTEPS 
RN  «  1.0 
HN  *  0.1*HN 
DTN  =  2 .0*HN 


IF (DTN  .LT. 

DTMIN)  H11  =  0.5*DTMIN 

f 

GO 

TO 

60 

50 

KGO  =  1 

i 

STEP  UNCHANGED  OR  INCREASED 

GO 

TO 

90 

60 

KGO  =  2 

i 

RESTART 

GO 

TO 

90 

70 

KGO  *  3 

1 

DT  <  DTMIN  (ERROR) 

GO 

TO 

90 

80 

KGO  =  4 

j 

STEP  DECREASE 

90 

RETURN 

END 

SUBROUTINE  ROTATE (N, INDEX) 

*** 

PURPOSE 

ROTATE  STACK  POINTERS 

*** 

INPUT 

N  NUMBER  OF  SUBVECTORS 

INDEX  ARRAY  OF  CURRENT  POINTERS 

*** 

OUTPUT 

INDEX  NEW  POINTERS 

INTEGER 

INDEX (1) 

IF (H-2) 30»10,20 

10 

I TEMP  = 

INDEX (2) 

INDEX (2) 

«  INDEX ( 1 ) 

INDEX (1) 
RETURN 

*  I TEMP 

20 

I TEMP  = 

INDEX (3) 

INDEX (3) 

=  INDEX (2) 

INDEX (2) 

•  INDEX ( 1) 

INDEX (1) 

=  I TEMP 

30 

RETURN 

END 

FUNCTION 

VIPDA (A»  B,  N) 

***  REAL  FUNCTION  VIPDA  ■  A  DOT  B,  A  AND  B  OF  LENGTH  N 
DOUBLE  PRECISION  ACCUMULATION 
RESULT  RETURNED  IN  SINGLE  PRECISION 


C 


DOUBLE  PRECISION  ACCUM 
REAL  A<1) ,B<1) 
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ACCUM  »  0 .DO 
DO  10  1=1, N 

ACCUM  =  ACCUM  +  DBLE(A(I) ) *DBLE(B(I) ) 
10  CONTINUE 

VI PDA  =  SNGL( ACCUM) 

RETURN 

END 


PARAMETERS  ASSOCIATED  WITH  CENTRAL  DIFFERENCE  INTEGRATOR 


COMMON  /CENTDF/ 

B  BTA,  1 

D  DECRS , DTMAX , DTMIN ,  I 

E  EPSDFQ, EPSHFQ, ERR, 

F  FIXSTP,  i 

H  HN, HNM1 , 

1  IDEC,  IDAMP, IDMPDG , IDTMX, IFORCE, INCRS, INCTRY , IRST,  ' 

2  11,12,13, 

K  KB AND,  5 

M  MAXDEG , MAX  STP , MAXVEC ,  S 

N  NCUTS, NDOUBL, NFACTS, NSTEP, NSTEPS , 

R  RN,  ! 

T  TIME,TMAX,  j 

U  UPNTR 

INTEGER  UPNTR ( 3)  j 

LOGICAL  DECRS, FIXSTP, IDAMP, IDMPDG, IFORCE, INCRS  j 


Appendix  B 

User  Written  Subroutines  with  Sample  Problem 

SUBROUTINE  DRIVER  IN  THIS  EXAMPLE  DRIVER  IS  THE  MAIN  PROGRAl 

—  PURPOSE  TO  PROVIDE  THE  DATA  NEEDED  TO  CALL  STINTC  (ELEMENT 

STINT/C END IF)  FOR  CENTRAL  DIFFERENCE  TIME  INTEGRATK 

—  INPUT  USER  SUPPLIED,  SEE  LIST  NEEDED  BELOW 


OUTPUT 


THE  ARRAYS 


LOGIC(20) 
INTGR(  20) 
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REALN(20) 
C ( I CORE) 
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WHERE 


LOGIC(l) 

s 

IGASP,  IF 

LOGIC (2) 

= 

I FORCE,  IF 

LOGIC(3) 

= 

IDISP,  IF 

LOGIC (4) 

IVEL,  IF  . 

LOGIC(5) 

= 

• 

FIXSTP,  IF 

LOGIC(6) 

= 

IDAMP,  IF 

LOGIC (7) 

s 

IDMPDG,  IF 

.TRUE.  DMGASP  I/O  USED 
.FALSE.  UNFORMATTED  FORTRAN  I/O 
.TRUE.  A  FORCING  FUNCTION  IS  USED 
.FALSE.  NO  FORCING  FUNCTION 
.TRUE.  INITIAL  DISPLACEMENTS 
.FALSE.  NO  INITIAL  DISPLACEMENTS 
TRUE.  INITIAL  VELOCITIES 
FALSE.  NO  INITIAL  VELOCITIES 
.TRUE.  USE  A  FIXED  TIME  STEP 
.FALSE.  VARIABLE  STEP  WILL  BE  USED 
.TRUE.  PROBLEM  HAS  DAMPING 
.FALSE.  NO  DAMPING 
.TRUE.  DAMPING  MATRIX  IS  DIAGONAL 
.FALSE.  DAMPING  MATRIX  NONDIAGONAL 


NOTE  IF  IDMPDG  =  .TRUE. ,  THEN  IDAMP  HAS  THE  FOLLOWING  MEANING 

1)  IF  IDAMP  ■  .TRUE.  USER  WILL  SUPPLY  DAMPING  MATRIX  AS 

A  VECTOR,  TO  BE  READ  IN  BY  STINT 

2)  IF  IDAMP  =  .FALSE.  DAMPING  IS  SUPPLIED  BY  A  USER  ROUTINE 

THAT  COMPUTES  D*UD  AS  IN  THE  NONDIAGONAL  DAMPING  CASE 


INTGR(l)  *  MAXDEG ,  THE  NUMBER  OF  DEGREES  OF  FREEDOM 
INTGR(2)  *  IUNIT,  THE  EXTERNAL  MASS  STORAGE  UNIT  NUMBER 
THAT  CONTAINS  THE  MASS  MATRIX,  DAMPING 
MATRIX  (IF  IDMPDG  AND  IDAMP  ARE  .TRUE.), 
INITIAL  DISPLACEMENT  AND  VELOCITY  (IF 
PRESENT) .  ALL  QUANTITIES  ARE  VECTORS  OF 
LENGTH  MAXDEG. 

I NTGR ( 3 )  ■  KBAND,  THE  WIDTH  OF  THE  CONNECTIVITY  BAND  OF 

THE  STIFFNESS  MATRIX  (OPERATOR)  OVER  WHICH 
THE  ERROR  IS  TO  BE  NORMALIZED. 

NOTE  KBAND  =  1  IMPLIES  EACH  DOF  IS  TREATED 
EQUALLY 

KBAND  =  MAXDEG/ 10  +  2,  RECOMMENDED 
(WHERE  MAXDEG/ 10  IS  INTEGER  ARITH¬ 
METIC) 


REALN(l)  = 
REALN(2)  = 

REALN(  3)  = 
REALNU)  « 


REALN(5)  « 


REALM (6)  « 


TIME,  THE  STARTING  VALUE  OF  TIME  FOR  THIS  RUN 
TMAX ,  THE  VALUE  OF  TIME  AT  WHICH  THE  INTEGRATOR 
WILL  STOP 

DTMIN,  THE  MINIMUM  TIME  STEP 
DTMAX,  THE  MAXIMUM  TIME  STEP 

(NOTE,  IF  FIXSTP  *  .TRUE.  STINTC  WILL  USE 
TIME  STEP  *  DTMIN  AND  SET  DTMAX'DTMIN) 
SAMPLES/CYCLE,  THE  NUMBER  OF  SAMPLES/CYCLE 
FOR  THE  DOMINANT  FREQUENCY  COMPONENT, 

MUST  BE  PI  OR  GREATER  FOR  STABILITY. 

10  TO  20  SAMPLES  IS  SUGGESTED 
4.0  IS  DEFAULT,  IF  REALN(5)  .LT.  PI 
SAMPLES/CYCLE,  THE  NUMBER  OF  SAMPLES/CYCLE 
FOR  THE  HIGHEST  APPARENT  FREQUENCY 
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COMPONENT,  MUST  BE  PI  OR  GREATER 
FOR  STABILITY. 

4PI  IS  INCREDIBILY  STABLE 
4.0  IS  DEFAULT,  IF  REALN( 6)  .LT.  PI 
REALN(7)  «  ALFA,  WHERE  THE  VALUE  OF  ALFA  LIES  BETWEEN  0.< 
AND  1.0.  SUGGEST  ALFA  NEAR  1.0  FOR  LIGHT  DAM1 
ING  AND  NEAR  0.20  FOR  CRITICAL  DAMPING  (SEE 
REFERENCES  3  &  4) .  ONLY  USED 
FOR  NONDIAGONAL  DAMPING 


ICORE  IS  THE  SIZE  OF  THE  DATA  ARRAY  C  THAT  STINTC  NEEDS 
TO  DO  ITS  THING.  THE  USER  MUST  INSURE  THAT  THIS 
ALLOCATION  IS  AVAILABLE. 


C 

c 


c 


c 


c 


c 


ICORE  -  (3 
WHERE 


+  3*L1  +  L2  +  L3  +  L4)*MAXDEG 


LI 

« 

3 

FOR 

VARIABLE 

L2 

■ 

1 

FOR 

DAMPING, 

L3 

m 

1 

FOR 

DIAGONAL 

L4 

a 

1 

FOR 

VARIABLE 

STEP,  1  FOR  FIXED  STEP 
0  OTHERWISE 
DAMPING,  0  OTHERWISE 
STEP,  0  OTHERWISE 


LOGICAL  LOGIC (20) 
INTEGER  INTGR( 20) 
REAL  REALN(20) 
COMMON  C(100) 


REAL  SM (2) , VZ ( 2 ) 


LOGIC (1)  -  .FALSE. 
LOGIC(2)  =  .TRUE. 
LOGIC (3)  »  .FALSE. 
LOGIC(4)  a  .TRUE. 
LOGIC (5)  »  .FALSE. 
LOGIC(6)  -  .TRUE. 
LOGIC (7)  »  .FALSE. 


I  UNFORMATTED  FORTRAN  I/O 
!  FORCING  FUNCTION 
!  NO  INITIAL  DISPLACEMENT 
l  INITIAL  VELOCITY 
!  VARIABLE  TIME  INCREMENT 
1  DAMPING  PRESENT 
!  NONDIAGONAL  DAMPING 


INTGR(l)  =  2 

! 

2  D.O.F. 

INTGR(2)  »  1 

1 

DATA  ON  UNIT 

1 

INTGR(3)  =  1 

! 

KBAND  *  i 

REALN(l)  a  0.0 

j 

START  TIME 

REALN( 2)  =  1.0 

i 

FINISH  TIME 

REALN( 3)  =  0.00001 

! 

MINIMUM  TIME 

INCREMENT 

REALM ( 4 )  -  1.0 

1 

MAXIMUM  TIME 

INCREMENT 

REALN( 5)  »  50.0 

1 

SAMPLE/CYCLE 

DOMINANT  FRE' 

REALN( 6)  =  4.0 

! 

SAMPLE/CYCLE 

HIGHEST  FREQ 

REALN(7)  -1.0 

1 

ALFA  -  DAMPING  PARAMETER 

VZ(1)  -  100.0 

1 

INITIAL  VELOCITY  D.O.F.  1 

VZ(2)  -  0.0 

I 

INITIAL  VELOCITY  D.O.F.  2 

SM ( 1 )  -  1.0 

! 

MASS  D.O.F. 

1 

SM(2)  -  1.0 

1 

MASS  D.O.F. 

2 

C 


OPEN ( UNIT-1 ,  FORM- ' UNFORMATTED ' ,  TYPE- ' SCRATCH ' ) 

WRITE  (1)  (SM(I) ,1-1,2)  1  UNFORMATTED  OUTPUT  TO  UNIT  1 
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WRITE  (1)  (VZ(I) ,1=1,2)  !  UNFORMATTED  OUTPUT  TO  UNIT  1 

REWIND  1 
C 

CALL  STINTC ( LOGIC, INTGR, REALM, C) 

C 

WRITE (6,1)  r 

1  FORMAT(////, '  ***  THAT1 'S  ALL  FOLKS  ***') 

END 


SUBROUTINE  FORCE (MAXDEG , TIME , F) 

-  PURPOSE  TO  PROVIDE  THE  USER  SUPPLIED  FORCING  FUNCTION,  F 

-  INPUT  MAXDEG, TIME 

-  OUTPUT  F 

-  DECLARATIONS  INTEGER  MAXDEG 

REAL  TIME, F ( 1) 

-  DEFINITIONS 

MAXDEG  *  NUMBER  OF  DEGREES  OF  FREED ON  (LENGTH  OF  F  ARRAY) 

TIME  =  THE  TIME  AT  WHICH  THE  FORCE  IS  TO  BE  CALCULATED 

F  =  THE  VECTOR  CONTAINING  THE  FORCE  AT  TIME.  F  IS  OF 
LENGTH  MAXDEG. 


-  EXAMPLE 

-  D.O.F.  2  SQUARE  FORCE  =  3000.  0.5<TIME<0.55 

-  FORCE  =0.0  OTHERWISE 

REAL  F ( 1) 

F( 1)  =  0.0 
F(2)  =  0.0 

IF (TIME  .LT.  0.5)  GO  TO  10 

IF (TIME  .GT.  0.55)  GO  TO  10 

F(2)  =  3000. 

10  RETURN 
END 


SUBROUTINE  SFORCE (MAXDEG , U, FS) 

PURPOSE  TO  PROVIDE  THE  USER  SUPPLIED  STIFFNESS  FORCES,  (FS) , 
I.E.  (FS)  ■  IK1 * (U) 


INPUT 


MAXDEG, U 
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-  OUTPUT 


DECLARATIONS  INTEGER  MAXDEG 
REAL  U(l) rFS(l) 

DEFINITIONS 

MAXDEG  «  NUMBER  OF  DEGREES-OF-FREEDOM  (LENGTH  OF  U  All 
ARRAYS) 

U  *  THE  VECTOR  CONTAINING  THE  DISPLACEMENTS 
FS  =  THE  VECTOR  CONTAINING  THE  STIFFNESS  FORCES 


-  EXAMPLE 

REAL  U(l)  /FSd) 

-  SPRING  RATE  COEFFICIENTS 

DATA  SK1/100Q./ / SKC/100 ./ f SK2/1000 ./ 

-  SPRING  (STIFFNESS)  FORCES  FOR  D.O.F.  1  AND  D.O.F.  2 

FS ( 1 )  *  SKl*TANH ( U ( 1 ) )  +  SKC* (U(l) -U(2) ) 

FS ( 2 )  -  SK2*SINH(U(2) )  +  SKC* ( U( 2) -U ( 1) ) 

RETURN 

END 


SUBROUTINE  DFORCE (MAXDEG f UD, FD) 

PURPOSE  TO  PROVIDE  THE  USER  SUPPLIED  DAMPING  FORCES,  (FD) , 
I.E.  (FD)  -  ID) * (UD) 


-  INPUT 


MAXDEG, UD 


-  OUTPUT 


DECLARATIONS  INTEGER  MAXDEG 

REAL  UD( 1) , FD ( 1 ) 

DEFINITIONS 

MAXDEG  »  NUMBER  OF  DEGREES-OF-FREEDOM  (LENGTH  OF  UD  AND 
ARRAYS) 

UD  «  THE  VECTOR  CONTAINING  THE.  VELOCITIES 
FD  -  THE  VECTOR  CONTAINING  THE  DAMPING  FORCES 


-  EXAMPLE 

REAL  UD ( 1 ) , FD ( 1 ) 
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DATA  D  /5.0/  !  DAMPING  VALUE 


FD  (  1  )  - 

DMUD(l)  -  UD(2)  ) 

!  D.O.F.  1 

DAMPING 

FORCE 

FD(2)  a 

D*(UD<2)  -  UD ( 1 ) ) 

!  D.O.F.  2 

DAMPING 

FORCE 

RETURN 

END 

< 

i 

SUBROUTINE  OUTPUT ( I PRT, MAXDEG , TIME , U, UD) 

-  PURPOSE  TO  SUPPLY  THE  DISPLACEMENT  AND  VELOCITY  AT  EACH  TIME 

STEP.  THE  USER  CAN  THEN  USE  THIS  DATA  FOR  PRINTING 
OR  PLOTTING  OUTPUT 

-  INPUT  I PRT, MAXDEG, TIME, U,UD 

-  OUTPUT  NONE 

-  DECLARATIONS  INTEGER  I PRT, MAXDEG 

REAL  TIME, U(l) , UD ( 1 ) 

-  DEFINITIONS 

I PRT  =  -2  STINT  HAS  ERRORED  OFF  NO  MORE  OUTPUT 
■  0  AT  ALL  TIMES  OTHER  THAN  BEGINNING  TIME 
»  2  AT  TIME  «  BEGINNING  TIME 
«  10  AT  THE  LAST  TIME  (PROBLEM  IS  FINISHED) 
NOTE,  THE  IPRT»2  FLAG  CAN  BE  USED  TO  INITIALIZE 
AND/OR  SET  UP  PRINTING  AND/OR  PLOTTING 
BUFFERS, I/O, ETC. 

MAXDEG  *  NUMBER  OF  DEGREES  OF  FREEDOM 
(LENGTH  OF  U  AND  UD  ARRAYS) 

TIME  *  THE  VALUE  OF  THE  CURRENT  TIME 

U  a  THE  CURRENT  DISPLACEMENT  VECTOR  (MAXDEG  VALUES) 

UD  *  THE  CURRENT  VELOCITY  VECTOR  (MAXDEG  VALUES) 


-  EXAMPLE 

REAL  U(l) ,UD(1) 

DATA  ICONT/O/ 

I F ( I PRT  .EQ.  -2)  GO  TO  30 

-  OUTPUT  RESPONSE  TO  UNIT  10  FOR  POST  PROCESSOR  PLOTS 

IF (I PRT  .EQ.  2)  OPEN(UNIT=10,  FORM- ' UNFORMATTED ' ,  TYPE* ' SCRATCH ' ) 
WRITE(IO)  TIME, (U(I) ,1-1, MAXDEG) , (UD(I> ,1-1, MAXDEG) 

-  PRINT  INITIAL  CONDITIONS 

I F ( I PRT  .EQ.  2)  WRITE(6 ,1 ) 

1  FORMAT (1H1) 

IF(IPRT  .EQ.  2)  GO  TO  10 

I F ( I PRT  .EQ.  10)  GO  TO  10 
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ICONT  *  ICONT  +  1 

C  -  PRINT -AT  EVERY  100-TH  TIME  INCREMENT 

I F ( MOD ( I CONT ,100)  .NE.  0)  GO  TO  20 

C 

10  WRITE (6, 2)  TIME, ( I , U( I ) , UD ( I ) , 1=1 , MAXDEG) 

2  FORMAT ( / , 

1  '  TIME  = ' , E10 . 3 ,/ , 

2  '  D.O.F.  DISPLACEMENT  VELOCITY',/, 

3  (7X,I1,6X,E10.3,6X,E10.3>) 

WRITE (6,3) 

3  FORMAT (/) 

C 

20  RETURN 

C  -  ERROR  EXIT 

30  WRITE (6,4) 

4  FORMAT ( // , '  *****  C END IF  HAS  ERRORED  OFF  ♦****•) 
STOP 

END 


Appendix  C 

Sample  Problem  Results 


WELCOME  TO  STINTC  THE  STAND-ALONE  CENTRAL  DIFFERENCE  TIME  INTEGRATE 
STINTC  REQUIRES  28  WORDS  OF  STORAGE 

6  DECEMBER  1979  VERSION 


TIME  *  0.000E+00 


D 

.  O.  F. 

DISPLACEMENT 

VELOCITY 

1 

0. 

000E+00 

0.100E+03 

2 

0. 

000E+00 

0.000E+00 

$$$ 

STEP 

SIZE 

INCREASED 

TO 

0.1500E-03 

AT 

0 . 1100E-02 

$$$ 

STEP 

SIZE 

INCREASED 

TO 

0. 2250 E- 03 

AT 

0.1850E-02 

$$$ 

STEP 

SIZE 

INCREASED 

TO 

0.3375E-03 

AT 

0.2975E-02 

$$$ 

STEP 

SIZE 

INCREASED 

TO 

0.5063E-03 

AT 

0.4663E-02 

$$$ 

STEP 

SIZE 

INCREASED 

TO 

0.7594E-03 

AT 

0.7194E-02 

$$$ 

STEP 

SIZE 

INCREASED 

TO 

0.1139E-02 

AT 

0 . 1099E-01 

$$$ 

STEP 

SIZE 

INCREASED 

TO 

0.1709E-02 

AT 

0.1669E-01 

$$$ 

STEP 

SIZE 

INCREASED 

TO 

0.2563E-02 

AT 

0.2523E-01 

sss 

STEP 

SIZE 

INCREASED 

TO 

0.3844E-02 

AT 

0 . 3804E-01 

$$$ 

STEP 

SIZE 

INCREASED 

TO 

0.5767E-02 

AT 

0.5727E-01 

$$$ 

STEP 

SIZE 

INCREASED 

TO 

0.8650E-02 

AT 

0.1149E+00 

TIME  »  0 . 409E+00 


CENTRAL  DIFFERENCE  TIME  INTEGRATOR 


D.O.F. 

1 

2 


DISPLACEMENT  VELOCITY 
0.739E+00  -0.5I6E+02 

-0.295E+00  -0.862E+01 


s  $$$ 

STEP 

$$$ 

STEP 

$$$ 

STEP 

$$$ 

STEP 

SIZE  DECREASED 
SIZE  INCREASED 
SIZE  DECREASED 
SIZE  INCREASED 


TO  0.6107E-02 
TO  0 . 8145E-02 
TO  0 . 5430E-02 
TO  0 .7701 E- 02 


AT  0.7723E+00 
AT  0.B273E+00 
AT  0 . 8517E+00 
AT  0.9169E+00 


TIME  =  0.100E+01 

D.O.F.  DISPLACEMENT  VELOCITY 

1  -0.199E-01  0.103E+02 

2  0.271E+00  0.312E+02 


AVERAGE  TIME  STEP  0.56818E-02 

NUMBER  OF  STEP  INCREASES  13 

NUMBER  OF  STEP  DECREASES  2 

NUMBER  OF  FACTORIZATIONS  0 

NUMBER  OF  TIME  STEPS  177 

NUMBER  OF  SOLUTIONS  179 


DT  OCCURRENCES  IN  THE  RANGES  INDICATED 


FROM 

10.0 

TO 

20.0 

TIMES 

DTMIN, 

DT 

OCCURRENCES 

WAS 

16 

FROM 

20.0 

TO 

30,0 

TIMES 

DTMIN, 

DT 

OCCURRENCES 

WAS 

5 

FROM 

30.0 

TO 

40.0 

TIMES 

DTMIN, 

DT 

OCCURRENCES 

WAS 

5 

FROM 

50.0 

TO 

60,0 

TIMES 

DTMIN, 

DT 

OCCURRENCES 

WAS 

5 

FROM 

70.0 

TO 

80.0 

TIMES 

DTMIN, 

DT 

OCCURRENCES 

WAS 

5 

FROM 

100.0 

TO 

200.0 

TIMES 

DTMIN, 

DT 

OCCURRENCES 

WAS 

10 

FROM 

200.0 

TO 

300.0 

TIMES 

DTMIN, 

DT 

OCCURRENCES 

WAS 

5 

FROM 

300.0 

TO 

400.0 

TIMES 

DTMIN, 

DT 

OCCURRENCES 

WAS 

5 

FROM 

500.0 

TO 

600.0 

TIMES 

DTMIN, 

DT 

OCCURRENCES 

WAS 

22 

FROM 

600.0 

TO 

700.0 

TIMES 

DTMIN, 

DT 

OCCURRENCES 

WAS 

10 

FROM 

700.0 

TO 

800.0 

TIMES 

DTMIN, 

DT 

OCCURRENCES 

WAS 

10 

FROM 

800.0 

TO 

900.0 

TIMES 

DTMIN, 

DT 

OCCURRENCES 

WAS 

79 
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***  THAT'S  ALL  FOLKS  *** 


L _ I 


Figure  2.  Time  Integrator  Functional  Outline 


INCREMENT  -  SECONDS 


DISPLACEMENT 


O 


Figure  5.  Sample  Problem  -  D.O.F.  1  Displacement  History 


I 

I 


LEGEND 


0.2  0.4  0.6  0.8 


TINE  -  SECONDS 


Sample  Problem  -  D.O.F.  1  Velocity  History 


TIME  -  SECONDS 


Figure  8.  Sample  Problem  -  D.O.F.  2  Velocity  History 


